Plaquette valence-bond ordering in J x — J 2 Heisenberg antiferromagnet on the honeycomb lattice 
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We study S — 1/2 Heisenberg model on the honeycomb lattice with first and second neighbor antiferromag- 
netic exchange ( J\ — Ji model), employing exact diagonalization in both S z = basis and nearest neighbor 
singlet valence bond (NNVB) basis. We find that for 0.2 < J2/J1 < 0.3, NNVB basis gives a proper de- 
scription of the ground state in comparison with the exact results. By analysing the dimer-dimer as well as 
plaquette-plaquette correlations and also defining appropriate structure factors, we investigate possible symme- 
try breaking states as the candidates for the ground state in the frustrated region. We provide body of evidences 
in favor of plaquette valence bond ordering for 0.2 < J2/J1 < 0.3. By further increasing the ratio J2/J1, this 
state undergoes a transition to the staggered dimerized state. 

PACS numbers: 75.10.Jm 75.10.Kt, 75.40.Mg 



I. INTRODUCTION 

Quantum spin liquid (QSL) is non magnetic state of a cor- 
related matter for which there is no broken symmetry in the 
spin part of the ground state wave function. Hence the lo- 
cal magnetic moments remain disordered down to absolute 
zero (T = 0)|Q]|. The quantum ground state for QSL can be 
expressed as the superposition of many different configura- 
tions, such as linear combinations of the short range singlet 
valence bonds. This state is called resonating valence bond 
(RVB), originally proposed by Fazekas and Anderson as the 
ground state of the Heisenberg model on the triangular lat- 
tice 02J]. The singlet bonds in RVB state can be considered 
as pre-formed Cooper pairs, which under suitable conditions ( 
i.e. hole doping) may coherently propagate throughout the 
system, hence give rise to superconductivity [3]. 

Many strongly correlated systems are well described by 
Hubbard Hamiltonian whose ground state for large on site 
coulomb interaction is the Mott insulating state. In this state 
the electrons are localized on the atoms, nevertheless local 
charge fluctuations induce an Anti-ferromagnetic (AF) ex- 
change interaction between the spins of the electrons. Hence, 
AF Heisenberg model is an effective Hamiltonian for describ- 
ing the low energy excitations of the Mott insulators [4]. It has 
been proved that in ID, the Ground state of the Heisenberg 
model is a gapless (critical) spin liquid for S = 1/2-chain, 
while it is gapped spin liquid for S = 1-chain [5]. Finding 
the realizations of QSL in two and three dimensions has been 
the subject of many researches in recent years [1]. Quantum 
fluctuations as well as frustration may destroy the long range 
magnetic order in spin systems. When a spin system is frus- 
trated, it can not find a spin configuration to fully satisfy the 
interaction between each pair of spins. There are two mech- 
anisms for frustration: (i) Geometrical frustration, where the 
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lattice geometry is such that, it is not possible to minimize the 
interaction energy of all bonds at the same time, e.g. in trian- 
gular or Kagome lattice in 2D and pyrochlore lattice in 3D Jg] . 
(ii) When there are several competing exchange interactions, 
such as competition between first and second neighboring AF 
exchange interactions ( J\ — J2 model). Since the quantum 
fluctuations are larger in 2D, many attempts to find QSL are 
focused on the quasi two-dimensional frustrated spin systems 
with S = 1/2 0]. 

Two dimensional Heisenberg antiferromagnets apart from 
their own importance [4], received intensive attention in the 
context of layered high-T c superconducting (HTSC) mate- 
rial yfl. The ground state of S = 1/2 Heisenberg model with 
AF nearest neighbor (NN) exchange coupling on 2D bipartite 
lattices has been shown to be Neel ordered J8 H14I1 . Addition 
of next nearest neighbor (NNN) AF interactions frustrates the 
system and gradually destroys the (Neel) order. Since AF ex- 
change interaction encourages the singlet formation, the quan- 
tum ground state of AF Heisenberg model can be expressed in 
term of over complete set of valence bond (VB) basis which 
represent a total spin singlet state B15I1 . 

The J\ — J2 AF Heisenberg model on square lattice has 
been extensively studied and various VB states have been pro- 
posed to describe its disorder regime [16]. One example of 
such quantum states is nearest neighbor RVB (NNRVB) rep- 
resenting a spin liquid, which breaks neither translational nor 
rotational symmetries. However, in highly frustrated regime 
where the ground state is classically disordered and the SU(2) 
symmetry of the Hamiltonian is restored, there is no theo- 
rem to prevent breaking of the lattice translational symme- 
try. Therefore, in spite of earlier proposal of states with no 
symmetry breaking B17[]. states which break translational sym- 
metry were proposed IU8H20I1 . One example is the staggered 
dimerized state which breaks translational and rotational sym- 
metries of the lattice. Another candidate, that has been pro- 
posed recently for some lattices, is plaquette RVB (PRVB) 
wave functions in which the resonance of VBs is limited to 
one plaquette J2ll - l25ll . PRVB state breaks the translational 
symmetry, while preserves the rotational symmetry of the lat- 
tice. 




FIG. 1: The bipartite honeycomb lattice. Two sublattices are marked 
by black and white circles. Nearest neighbor lattice points are con- 
nected with solid lines and next to nearest neighbor lattice points are 
connected with dashed lines. Red arrows show the two primitive lat- 
tice vectors. 



Recent fabrication of graphene monolayer and also mag- 
netic compounds with quasi 2D honeycomb structure, has 
brought the honeycomb lattice to attention of the physicists 
from both experimental and theoretical points of view. Hon- 
eycomb lattice does have coordination number equal to three, 
which is minimum among two-dimensional lattices. In the 
case of Heisenberg model on honeycomb lattice, the small 
number of neighboring interactions enhances the quantum 
fluctuations and therefore seems to be promising system to ex- 
plore spin liquid states. Honeycomb lattice is a bipartite lattice 
composed of two interlacing triangular sublattices (Fig. [T]). 
The unit cell of this non-bravais lattice contains two sites and 
lattice is constructed by two lattice vectors of the triangular 
bravais lattice. The non-bravais character of lattice results in 
more exotic aspects that can not be seen in square lattice or 
the other bravais lattices [26]. 

As some realizations of Heisenberg magnets on the hon- 
eycomb lattice, one can name recently discovered com- 
pounds such as InCu 2 / 3 V 1 /303 JH] and Na 3 Cu 2 Sb0 6 JH 
in which the Cu +2 ions in the copper-oxide layers form a 
two-dimensional S = 1/2 Heisenberg antiferromagnet on a 
honeycomb lattice, Bi 3 Mn 4 0i2(N03) (BMNO) in which the 
Mn +4 ions with S = 3/2 reside on the lattice points of weakly 
coupled honeycomb layers H29I1 . Replacing Mn +4 with V +4 
in this compound may realize the S = 1/2 Heisenberg model 
on honeycomb lattice. Also the recent progress in the field of 
ultracold atoms and trapping techniques [30] along with the 
ability to tune the interaction parameters via the Feshbach res- 
onance [31] can be thought of another way to realize Heisen- 
berg spins (of localized fermions) on a honeycomb optical lat- 
tice. 

Two recent achievements has raised the hope of finding 
QSL in honeycomb geometries. One, is the large scale quan- 
tum Monte Carlo simulation of the half-filled Hubbard model 
on the honeycomb lattice, which results a spin liquid phase 
with finite spin gap for moderate values of on-site coulomb 
interaction (3.5 < U/t < 4.3). This phase is located be- 



tween the semi-metallic phase characterized by massless dirac 
fermions (U/t < 3.4) and the AF-Mott insulating phase for 
U/t > 4.3 13211 . The other is the experimentally observed 
spin liquid behavior in BMNO which remains magnetically 
disordered down to T — 0.4K, in spite of its high Curi-Weiss 
temperature T cw ~ -257K J29J1 . 

Motivated by the above considerations, in this paper we in- 
vestigate the ground state properties of J\ — J 2 Heisenberg 
model which is proposed for explaining the spin liquid behav- 
ior of BMNO. The paper is organized as follows. In section. ITT1 
we introduce the spin Hamiltonian and using diagonalization 
in nearest neighbor VB basis, we find evidence for spin liquid 
phase for a range of coupling constants. In section. iDTl we em- 
ploy the exact diagonalization in full Hilbert space of S z — 0. 
With exact wave-functions obtained in this manner, we cal- 
culate the dimer-dimer correlation function. We find two dif- 
ferent quantum phases in frustrate regime. In section. [TV] by 
introducing suitable structure factors, quantum phase transi- 
tion point is determined. At the end, in section. [V] we calcu- 
late plaquette-plaquette correlations which points to a possible 
PRVB state. Section. [VIlis devoted to discussion and conclu- 



II. MODEL HAMILTONIAN AND ITS GROUND STATE 
CANDIDATES 

3\ — J 2 AF Heisenberg Hamiltonian is defined by, 
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in which J\ > and J 2 > are AF exchange interactions 
between first and second neighboring spins, respectively. The 
first sum is limited to NN sites, while the second sum runs 
over the NNN lattice sites. Since the square lattice is con- 
nected with high T c superconducting materials, the studies of 
frustrated phases of spin models has been usually limited to 
this lattice. Recently discovered magnetic materials with un- 
derlying honeycomb geometry is our motivation to study the 
above model on the honeycomb lattice. 

Using effective action approach to the frustrated Heisen- 
berg antiferromagnet in two dimensional system, Einarsson 
et. al. proved that the S = h disorder ground state on the hon- 
eycomb lattice is three-fold degenerate [33]. Fouetef. al. 113411 
provided some evidence for staggered dimerized (SD) state 
for S — 1/2 at J 2 / Ji = 0.4. Such state breaks the rotational 
lattice symmetry (C 3 ), while preserves its translational sym- 
metry (Fig. |2). They also speculated that spin liquid phase 
and PRVB phases can be stabilized for some range of J2/J1 
(Fig [3}. PRVB phase breaks the translational symmetry, but 
preserves the rotational symmetry of lattice. Alternatively, 
Read and Sachdev [35] used large N expansion method, pro- 
posed another ground state wave function which breaks both 
the translational and the rotational symmetries of the lattice 

(Figffll. 

In the classical limit (large spins), it has been shown that 
the ground state of the above model is Neel ordered for 





FIG. 2: Three degeneracy of staggered dimerized wave function on 
honeycomb lattice. 



FIG. 4: Three degeneracy of wave function proposed by Read and 
Sachdev on honeycomb lattice. 




FIG. 3: Three degeneracy of Plaquette Valence Bond wave function 
on honeycomb lattice. 



J2/J1 < 1/6 while for J2/J1 > 1/6 the ground state con- 
sists of an infinitely degenerate set of spiral states character- 
ized by spiral wave vectors q [36]. Okumura et al [37] used a 
combination of low temperature expansion and Monte Carlo 
simulation showed that thermal fluctuations can lift the huge 
degeneracy of the ground state, leading to a state with bro- 
ken C3 symmetry of honeycomb lattice. According to their 
finding, in the vicinity of Neel phase boundary, the energy 
scale associated with the thermal order-by-disorder becomes 
so small that exotic spin liquid behavior, such as ring-liquid 
or pancake-liquid can emerge. Mulder et. al. argued that 
taking the quantum fluctuations into account, some specific 
wave vectors in this manifold are picked as the ground state 
- a manifestation of order by disorder mechanism. They find 
for S = 1/2 quantum fluctuations are strong enough to de- 
stroy the spiral order and stabilize the valence bond solid with 
staggered ordering [38]. Our aim in this paper is to study the 
ground state of model (Q]i using exact diagonalization in both 
£2 and nearest neighbor valence bond (NNVB) basis. 



III. DIAGONALIZATION IN NNVB BASIS 

The valence bond states are a subset of S Z =Q basis with 
total spin magnitude S 2 equal to zero. In this section we show 



that the ground state of the J\ — J2 Heisenberg model in the 

frustrated regime, where there is no long range order, can be 

very well approximated in terms of states in NNVB subspace. 

Let us expand the ground state wave function in terms of 

NNVB states as 



H'q) = ^2a(c a )\c a ), 



(2) 



where \c a ) denotes all possible configurations a of NNVBs: 
l c a> = II (Hh-hM- ( 3 ) 

(i,j)<£a 

First, we have to enumerate the basis \c a ) to construct a nu- 
merical representation of the Hamiltonian matrix in this ba- 
sis. To determine the basis, the exact Pfaffian representation 
of the RVB wave function is employed [39]. In this method 
one expresses the RVB wave function as the Pfaffian of an an- 
tisymmetric matrix whose dimension is equal to the number 
of lattice points. The NNVB basis is much smaller than the 
whole S z = basis, so that the Hamiltonian matrix can be 
fully diagonalized with standard library routines. Note that 
since the NNVB states (|c Q )) are not orthonormal, one needs 
to solve the generalized eigen-value problem 

det[H - EO] = 0, 

where O — {cp\c a ) denotes the overlap matrix of the NNVB 
configurations. 

In the upper panel of Fig. [5] we have compared ground state 
energies obtained in the NNVB basis, and those obtained by 
numerically exact diagonalization in the S z = basis ver- 
sus J%l J\. In the middle panel we show the relative error in 
the ground state energy and the lower panel shows the over- 
lap of the exact ground state wavefunction with the ground 
state obtained within NNVB basis set. As can be seen in 
this figure, the agreement between the two sets of energies 
for J2/J1 € ]0.2, 0.3[ is remarkable. Since the NNVB ba- 
sis is not complete, the large error obtained by NNVB ba- 
sis for J2/J1 < 0.2 and J2/J1 > 0.3 can be attributed to 
the fact that longer range valence bonds start to contribute. 
For J2/J1 < 0.2, where there is Neel order in the ground 
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FIG. 5: Up: The comparison between the exact ground state energy 
evaluated using exact diaganolazition in ^=0 basis (squares) and di- 
agonlization in NNVB basis (circles) as a function of J2/J1 ■ middle: 
The relative errors between the ground state energies obtained by the 
two basis sets denned as (E^ VB -E$ I> )/(E$ U ). Down: The over- 
lap of the exact GS wavefunction lipo } with the GS wavefunction 
obtained in NNVB basis |V>o™ VB )- 



FIG. 6: Probability distribution of the relative amplitudes of VB's in 
the ground state wave-function for J2/J1 = 0.2, 0.4 and cluster size 

N = 54. 



plitude. Thus it can be represented as: 



\HVB) = A^ i \c a ), 



(4) 



in which A = 



is the normalization coef- 



ficient. In order to get preliminary insight into the nature of 
the ground state obtained in NNVB basis, we compare the 
distribution of the amplitudes of VB configurations in eq. (0 
to the uniform amplitude A of the above RVB liquid state. 
For this purpose we define the ratio w(c a ) — 'v , and look 
at the distribution of relative weights, w. If the ground state 
has the characteristics of a RVB spin liquid, this distribution 
is expected to be sharply peaked around w = 1. In Fig. [6] 
we plot the probability distribution function (PDF) of the rel- 
ative amplitudes w for J2/J1 = 0.2 and J2/J1 = 0.4 in a 
cluster of N = 54 spins. As can be seen in this figure, for 
J2/J1 = 0.2 the PDF is narrower relative to J2/J1 = 0.4, 
which implies that the ground state for J2/J1 — 0.2 is more 
similar to spin liquid state. The broader distribution of ampli- 
tudes for J2/J1 > 0.4 on the other hand, indicates significant 
deviation from QSL behavior, which can be considered as a 
sign of symmetry breaking. To investigate the effect of finite 
lattice geometry on the ground state, we compare PDFs of rel- 



state, it was shown that long-ranged VB states have remark- 
able contribution in the ground state wave function [40]. Start- 
ing from J2/J1 = 0, the Neel order is destroyed by increas- 
ing frustration strength up to J 2/ J\ ~ 0.2. At this point, the 
spin-spin correlations will become short ranged and the na- 
ture of the ground state can be accurately captured by NNVB 
wave functions. For J2/J1 > 0.3, the frustrating second 
neighbor AF exchange coupling J2 induces singlet formation 
between next nearest neighbors (NNNVB's) which compete 
with NNVB's. Therefore in this region the NNVB basis is in- 
sufficient to capture the true ground state. Furthermore, upon 
increasing J2/J1 beyond 0.35, as will be shown shortly, the 
nearest neighbor singlets in the VB states (dimers) will start 
to become correlated. 

The RVB wave function (QSL) is defined as linear super- 
position of all possible VB configurations with the same am- 
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FIG. 7: (Color online) PDF of relative amplitude w for J2/J1 = 
0.2, 0.35 for (a) N = 50 and (b) N = 54. (c) Width and (d) skew- 
ness of PDF vs. J2/J1 for N = 50, 54. 
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TABLE I: {P a P a ') — {Pa} 2 for a fixed and a' = /3,7, S (Fig. |8). Three indices 1, 2, 3 refer to three degenerate states (cf. Figs. [2] [3] |4) 
which become orthogonal to each other in the thermodynamic limit for pure SD and RS states. Subscript avg denotes average over these 
three possible degeneracies. Since the three degenerate PL states are not orthogonal to each other in the thermodynamic limit, in this case the 
correlations have been extracted from finite size scaling of the numeric data. Digits in the parenthesis denote errors in the last digit. 



ative amplitudes for two cluster sizes N — 50 and N = 54, 
and J2/J1 — 0.2 and J2/J1 = 0.35. These are shown in 
the two top panels of Fig. Q It can be easily seen that for 
the lattice size N = 50, the PDFs are more symmetric than 
N = 54. This can be attributed to the fact that, in contrast to 
the cluster with N = 50 lattice points, the plaquette wave 
function can be fitted to the cluster with N = 54 subject 
to the periodic boundary condition. Therefore the symme- 
try breaking toward plaquette formation is more pronounced 
for N = 54. This signals a tendency for plaquette formation, 
provided it is compatible with the lattice geometry. Moreover 
for N — 54, increasing the value of J2 leads to emergence of 
a second peak at the right tail of PDF. The average of relative 
amplitudes of VB configurations taking part in the plaquette 
wave function turns out to be 1.56. Interestingly, position of 
the second peak is quite close to this values. Therefore, the 
second peak can be attributed to the amplification of plaque- 
tte character in the ground state, as a result of increasing the 
second neighbor exchange interaction. To quantify the vari- 
ations of PDF (p(w)) versus J2/J1, we compute its width 
given by the standard deviation (a = {(w — (ui)) 2 ) 1 / 2 ) and 
also its skewness defined by u^-wj ) _ p ane i ( c ) f pjg [7] 
is plot of width versus J2/J1 for N = 54, showing that the 
width of PDF increases monotonically from J2/J1 = 0.2 to 
J2/J1 — 0.4. The panel (d) in this figure, represents the rise 
of PDF skewness for N = 54 in terms of J2/J1, indicating 
that by increasing J2 the distribution functions for this size get 
more and more asymmetric, while for A^ = 50 skewness does 
not change remarkably. This suggests that the skewness can 
be considered as a heuristic indication of possible symmetry 
breaking. In summary, since the widths of PDFs considered 
here are finite for the interval 0.2 < J2/J1 < 0.35, the above 
arguments suggest that the ground state of J\ — J2 model in 
this interval is not a perfect QSL state. A more precise under- 
standing of the ground state properties, requires the study the 
correlation between dimers. This will be done in the following 
section. 



dimers for 0.2 < J2/J1 < 0.5. The dimer-dimer correlation 
is defined by, 

C(a,a') = 4(((S i .S i )(S fc .S,)) - ((S^S,)) 2 ), (5) 

where a' = (k, I), and a — (i,j) is the reference bond relative 
to which the correlations are calculated. Define the permuta- 
tion operator by, 



P ki =2(S fc .S,) + i, 



(6) 



in terms of which Eq. (O can be alternatively expressed as 
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In table U quantities C(a,a r ) for fixed a and a' = /3,7, S 
(Fig. [8]l are shown for three trial wave functions V'sd, V'pl 
and i/'RS, where RS, SD and PL stands for Read-Sachdev, 
staggered dimerized, and plaquette states, respectively. The 
expectation values of the operator (P a i ) for each of the three 
degenerate SD and RS states is —1 if bond a' is occupied 
by a dimer, and +1/2, otherwise. (P a ')a,vg in table [I] is ob- 
tained by averaging over three degenerate states correspond- 
ing to SD and RS trial wavefunction. For PL trial state, since 
the three degenerate configurations are not orthogonal in ther- 
modynamic limit, the averaging is not valid and we compute 
C(a, a') numerically by finite size scaling method. Fig. [9] 
gives a graphical representation of the correlations obtained in 
this way. Red and blue links denote positive and negative cor- 
relations, and their thickness is proportional to the magnitude 
of correlations with the reference bond a of Fig. [8] The refer- 
ence bond a is denoted by a double line in Fig. [9] It can be 



IV. EXACT DIMER-DIMER CORRELATIONS 



In this section, employing exact diagonalization we obtain 
the ground state in S z = basis. Using the exact wave 
function of the ground state, we calculate correlation between 



FIG. 8: The reference bond a, and three independent bonds f3, 7, 5. 




FIG. 9: Snapshots of correlations corresponding from left to right to 
(a) RS, (b) PL, and (c) SD trial states. Red and blue links correspond 
to positive and negative correlations with respect to reference bond 
a which has been denoted by double line. 




FIG. 10: The dimer-dimer correlation for honeycomb lattice with 
periodic boundary condition at J2 = 0.3. Red (Blue) lines denote 
positive (negative) correlation. The thickness of lines is proportional 
to the magnitude of correlations. The system size is N — 32. 



seen in this figure that the pattern for the sign of correlations 
in PL and RS states are identical, except for their magnitudes. 
Specifically, in the PL state, the magnitudes of correlations 
between the dimers which belong to the same hexagon as the 
reference bond, are stronger than the rest of dimers. But in the 
RS state all dimers have identical correlation with respect to 
the reference bond. 

In Fig. [10] and Fig. [TTJ we have shown the exact diago- 
nalization results for the dimer-dimer correlation functions at 
J2/J1 — 0.3 and J2/J1 — 0.4, respectively for a lattice with 
N = 32 sites, subject to periodic boundary conditions. Note 
that in order to implement symmetries of infinite lattice on fi- 
nite size systems, the size N is limited to specific numbers 
N = 24, 32, 42, 50, 54, . . .. Since the dimension of Hilbert 
space grows exponentially with N, the exact diagonalization 
in the whole S z = subspace is not feasible, and hence for 
N > 32 we carried out the calculations in NNVB basis. The 
correlations are computed with respect to the middle-bond in- 
dicated by double lines. Red bonds denote positive correla- 




FIG 11: The dimer-dimer correlation for honeycomb lattice with 
periodic boundary condition at J2 = 0.4. Red (Blue) lines denote 
positive (negative) correlation. The thickness of lines is proportional 
to the magnitude of correlations. The system size is TV = 32. 



tions, while the blue ones indicate negative correlations. The 
thickness of bonds are proportional to the magnitude of corre- 
lations. As can be seen in Fig. [TTJ for J2/J1 = 0.3, the corre- 
lations are decaying with distance (measured with respect to 
central bond). Comparing this correlation map with Fig.[9jb), 
shows remarkable similarity to the snapshot corresponding to 
PL ordering. 

The dimer-dimer correlation pattern at J2/J1 = 0.4 shown 
in Fig. [TTJ is entirely different. First the dimer-dimer corre- 
lations do not appreciably decay over the maximum distance 
displayed in the figure. Second, one can easily distinguish 
a staggered ordering pattern by noting that, correlations be- 
tween bonds parallel to the reference dimer are positive, while 
others are negatively correlated with the reference dimer. 
Comparing with Fig.[9](c), this correlation snapshot obviously 
suggests a staggered dimerized state at J2/J1 = 0.4. This is 
in agreement with previous study of Fouet and coworkers [ 34] . 
However, Fouet et. al. speculated that at J2/J1 = 0.3 the cor- 
relation pattern resembles a RS state. Based on qualitative 
symmetry consideration of short-range dimer-dimer correla- 
tions, Fouet and coworkers proposed the possibility of crystal 
of hexagon plaquettes as a candicate for the ground state in the 
0.3 < J2/J1 < 0.35 range [34]. In the following we demon- 
strate that for 0.2 < J2/J1 < 0.3, the dominant correlations 
are of PL type, rather than RS. 

For quantitative characterization of the nature of VB crys- 
talline state, we define the following structure factor: 



S\ = 






ex(a') C(a,a'), 



(8) 



where C(a, a') is given by Eq. (TTJl and £\(a') is the phase fac- 
tor, appropriately defined for each of the three states A =SD, 
PL, RS [24]. The phase factors Ssd,£pl and £rs are shown 
in Fig. [T_2] Since the sign of dimer-dimer correlations for PL 
and RS states are the same, their phase factors must be equal. 
Scaling behavior of S\ for a lattice with N sites and Nf, bonds 
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FIG. 12: Red and blue links correspond to +1, —1 phase factors, 
respectively, while the dashed links stand for 0. Reference dimer is 
identified by double link, (a) Phase factor convention for PL and RS 
states, (b) Phase factors for staggered dimerized state. 



is given by, 



N~ b ~ Cx + N- 



(9) 



Using the above phase factors, and the correlations C(a, a') 
given in table U we have calculated the corresponding C£° in 
table|II]for each of the three trial states ip\. 

In Fig.[l3]we have shown S*sd and 5pl versus J2/J1 for 
N = 24,42,54. Note that for N = 24, we have done ex- 
act diagonalization in S z = basis, while for N — 42, 54, 
NNVB basis has been used. Note that the NNVB calculations 
is valid only in the region 0.2 < J2/J1 < 0.35. With these 
three sizes, we have performed finite size scaling according 
to Eq. © to obtain Cg^ and C^ L . Fig. [Ola) and (b) show 
SD and PL structure factors, respectively, for various sizes. 
For values of J2/J1 > 0.35 we do not have reliable data for 
N = 42, 54. Hence we have not reported finite size scaling 
for these values of J2/J1. For 0.2 < J 2 /J\ < 0.35, the SD 
structure factor in (a) remains much smaller than C^> = 1/4 
(table [Il]i for all three sizes. On the other hand, the PL struc- 
ture factor in (b) can be extrapolated to a finite value between 
0.07 and 0.1 (empty squares) which are comparable with the 
exact value 0.125 of the pure PL state (tablefTTb. 

A sudden jump observed in Fig.[l3](a) and (b) suggests the 
existence of a first order phase transition from PL to SD state 
as one increases J2/J1 beyond a certain value between 0.35 
and 0.4. In panels (b) of Fig. [I3]the average ratio of structure 
factor (averaged over the range 0.2 < J2/J1 < 0.35) to the 
corresponding Cp\ is given by (Spl/N^/C^ « 0.66, while 
this ratio for RS state in the same region is (Sp^s / Nt,) / Cp^ s « 
0.08/0.375 sw 0.21. Hence for 0.2 < J 2 /Ji < 0.35 we 
expect the ground state to be dominated by plaquette valence 
bond order. 



TABLE II: Intensive structure factor in thermodynamic limit for the 
three trial states. 
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V'RS 
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FIG. 13: The structure factor computed for lattice with N = 24 (di- 
amond), TV = 42 (triangle) and N — 54 (square). Nb stand for 
the number of dimers. Structure factors correspond to (a) Staggered 
dimerized state, (b) Plaquette state. In (b), empty square with er- 
ror bar indicates extrapolation to infinite lattice size, (c) shows a 
finite size scaling according to Eq. ij9) for PL structure factors and 
J2/J1 = 0.3. As can be seen for a specific value of J2/J1 between 
0.35 and 0.4, there is a suddent incrase for SD structure factor, which 
is accompanied by a sudden decrease in the PL structure factor. 



In addition we calculated the exact value of (P a ') as a 
function of J2/J1 for N = 24 sites. For J2/J1 from 0.2 
to 0.35, the expectation value (P a ') increases monotonically 
from —0.21 to —0.12 and shows a sudden jumps to 0.001 for 
J2/J1 = 0.4. In view of the ra —0.1 value for the expec- 
tation value of permutation operator in the PL state (table |I|, 
the negative values in the range 0.2 < J2/J1 < 0.35 can be 
considered as an extra support in favor of plaquette valence 
bond solid in this regime. Guided by the above evidences for 
plaquette ordering in region 0.2 < J2/J1 < 0.35, in the next 
section we calculate the plaquette-plaquette correlation using 
exact diagonalization method. 




TABLE III: Relative plaquette-plaquette correlations for hexagons 
with different distance to reference hexagon. 



FIG. 14: Red (blue) circles denote positive (negative) correlations. 
The radius of circles are proportional to the value of plaquette- 
plaquette correlation. The reference plaquette is depicted by black 
dashed line, (a) The plaquette-plaquette correlation map calcu- 
lated for PL trial state. The exact plaquette-plaquette correlation 
evaluated using exact diagonalization on honeycomb lattice for (b) 
J2/J1 = 0.3 and (c) J2/J1 — 0.5. (d) PL-PL correlation pattern for 
iV = 54and J2/J1 =0.3. 



C(i,0)/C(0,0) 


1/>PL 


VWb 


V'ED 


i = l 


-0.12 


-0.04 


-0.08 


i = 2 


0.28 


0.13 


0.175 


i = 3 


-0.12 


-0.04 


-0.07 


i = 4 


-0.12 


-0.008 


-0.05 


i = 5 


0.28 


0.05 


0.11 



in S z = basis for TV = 24 sites at J 2 /J\ = 0.3,0.5, re- 
spectively. Comparison of panels (b) and (c) of Fig. [14] with 
panel (a), indicates substantial PL ordering at J2/J1 = 0.3. It 
is remarkable to note that even the ratio of strengths of pos- 
itive and negative correlations in (b) (~ 0.29 : 0.14) and (a) 
(0.35 : 0.13) agree with each other. This ratio in (c) becomes 
0.04 : 0.04 which significantly deviates from the correspond- 
ing value for a pure PL state (a). Fig. |T4fd) is the same as 
(b), for larger size N — 54 and J2/J1 = 0.3. As can be seen, 
for N = 54 too, a substantial plaquette correlation pattern 
can be observed. Moreover, for this size we calculated the 
PL-PL correlations for all hexagons (denoted by 1,2,3,4,5 in 
FigfPftd)) with respect to the reference one to which, number 
is assigned. In last column of table [III] we have listed the 
relative correlations, defined by C(i, 0)/C(0, 0), obtained by 
exact diagonalization in NNVB basis and compare them with 
corresponding values for PL and RVB states given in first and 
second columns, respectively. As can be seen, for RVB state 
this quantity decays rapidly with distance from the reference 
hexagon, while the exact data does not show such a rapid de- 
caying. This again is an evidence in favor of the PL nature of 
the ground state. 



V. PLAQUETTE ORDER IN FRUSTRATED REGIME 

A more direct tool to detect plaquette ordering in frustrated 
regime is to investigate the plaquette-plaquette correlation de- 
fined by, 



C( P ,q) = (Q P Q q )~(Q P ) 2 , 



(10) 



where p, q stand for different plaquette and IT (II -1 ) is the 
cyclic exchange operator which permutes six spins around a 
hexagon in clockwise (counter-clockwise) direction. This cor- 
relation function was introduced recently and has been used to 
investigate plaquette ordering in frustrated Heisenberg mag- 
nets on the checkerboard and square lattice ll2lT423l 12511 . 

In Fig. [T4ka) we have depicted the plaquette correlation 
in PL state. Red (blue) circles indicate positive (negative) 
correlations, with the radii of circles are proportional to the 
magnitude of correlation. The reference plaquette is marked 
with a dashed circle. Fig. [T4l-(b) and (c) represent the pla- 
quette correlation function obtained by exact diagonalization 



VI. CONCLUSION 

In summary, diagonalization of J\ — J2 antiferromagnet 
Heisenberg Hamiltonian on honeycomb lattice in both S z =0 
basis and NNVB basis show a striking agreement between 
these two approaches in the parameter range 0.2 < J2/J1 < 
0.3. Therefore, in this region the ground state can be well 
described in terms of the singlet bonds between the nearest 
neighbor spins. Analysis of the exact dimer-dimer correla- 
tions, structure factors, and also plaquette-plaquette correla- 
tions, suggests the existence of plaquette valence bond crystal 
in this range of couplings. In fact the emergence of such 
PL ordering can be attributed to the quantum fluctuations due 
to the tendency of second neighbors to form singlets. This 
study also reveals a phase transition from the plaquette or- 
dered to the staggered dimerized state at a point in the interval 
J2/J1 €]0.35, 0.4]. Similar results, regarding plaquette or- 
dering on the square lattice, have been previously obtained 
for J\ — J2 — J3 Heisenberg antiferromagnet in its maximally 
frustrated region, J2 + J3 ~ Ji/2, and for J2 < J3 12411 . Our 
results are in contrast with those obtained by QMC simulation 
of the Hubbard model in intermediate interaction regime [32], 
in a sense that QMC results in, a RVB liquid phase with no 



9 



broken symmetry for this region. Meng, et. al. have only cal- 
culated short range dimer-dimer correlation. Therefore this 
difference might be due to the fact that the geometry of finite 
clusters used in their simulation is not compatible with PL or- 
dering. Based on present study, we believe that it is necessary 
to consider lattice geometries compatible with PL state and to 
calculate the plaquette-plaquette correlation for precise deter- 
mination of the broken symmetries in the ground state. 
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